function [meanIsp,stdIsp,Isps,variousIsps] = calcIsps(filenames,fuelmass,plotThings)
% calc Isps for given filenames and return mean, standard deviation and
% listing of all the isps
% inputs: filenames (cell array of strings)
% fuelmass: mass of fuel in kg, if not specified will use 1kg 

if nargin <2
    fuelmass = 1;
    plotThings = 0;
elseif nargin < 3
    plotThings = 0;
    
end
g = 9.81;
twater = 0.16; % seconds to end of water thrust - estimate, used for one 
% method of removing influence of bottle pressurization on Isp

nfiles = length(filenames);

impulse = zeros(1,nfiles);
Isps = zeros(1,nfiles);

wh = waitbar(0,'Calculating Isps');
for i = 1:nfiles
%     if ~rem(i,2)
%         waitbar(i/nfiles,wh,'Calculating Isps');
%     end
    
    [t,T,nomForces] = loadThrustData(filenames{i},plotThings);
    
    % get rid of last two points, that set the end of teh data set to zero
    t = t(1:end-2); T = T(1:end-2);
    
    % build vectors for subtracting off
    nTs = length(t);
    halfForce = mean(nomForces);
    tendInd = find(t > twater,1,'first');
    linBase = linspace(nomForces(1),nomForces(2),nTs);
    stepLinBase = [linspace(nomForces(1),halfForce,tendInd),...
        linspace(halfForce,nomForces(2),nTs - tendInd)];
    
    Tinit = T - nomForces(1);
    Tfin = T - nomForces(2);
    
    Tlin = T - linBase';
    Tstep = T - stepLinBase';
    
    %zero anything that ends up being negative
    Tinit = Tinit.*(Tinit >= 0);
    Tfin = Tfin.*(Tfin >=0);
    Tlin = Tlin.*(Tlin >=0);
    Tstep = Tstep.*(Tstep >=0);
    
    impulse(i) = trapz(t,Tfin);
    
    Isps(i) = impulse(i)/(fuelmass*g);
    
    IspsBase(i) = trapz(t,T)/(fuelmass*g);
    IspsFin(i) = Isps(i);
    IspsInit(i) = trapz(t,Tinit)/(fuelmass*g);
    IspsLin(i) = trapz(t,Tlin)/(fuelmass*g);
    IspsStepLin(i) = trapz(t,Tstep)/(fuelmass*g);
end
close(wh);

meanIsp = mean(Isps);
stdIsp = std(Isps);

variousIsps.base = IspsBase;
variousIsps.fin = IspsFin;
variousIsps.init = IspsInit;
variousIsps.lin = IspsLin;
variousIsps.step = IspsStepLin;


Isps = IspsLin;  %actually output linear offset correction